%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Performs ICA decomposition on MEG data read into FieldTrip. %
% Last modified: July 18, 2014                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2013-2014, Michael J. Cheung
%
% This file is a part of the MEG & PLS Pipeline (MEGPLS). For more 
% details, see the documentation included with the software package.
%
% MEGPLS is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License version 2 as published by 
% the Free Software Foundation. This program is distributed in the hope 
% that it will be useful, but WITHOUT ANY WARRANTY; without even the 
% implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along
% with this program. If not, you can download the license here: 
% <http://www.gnu.org/licenses/old-licenses/gpl-2.0>.


function MEGpipeline_RunICA(FTcfg, InputParams, RunChoice)


% Make sure java paths added:
UseProgBar = CheckParforJavaPaths;


% Load InputParams:
name  = InputParams.name;
paths = InputParams.paths;


% Set up errorlog and diary:
if exist('ErrorLog_ICA.txt', 'file')
    system('rm ErrorLog_ICA.txt');
end
if exist('Diary_ICA.txt', 'file')
    system('rm Diary_ICA.txt');
end

diary Diary_ICA.txt
ErrLog = fopen('ErrorLog_ICA.txt', 'a');



%==========================%
% RESAMPLE DATA & RUN ICA: %
%==========================%


for g = 1:length(name.GroupID)
    NumSubj = length(name.SubjID{g});
    
    if UseProgBar == 1
        ppm = ParforProgMon...
            (['RUNNING ICA DECOMPOSITION:  ',name.GroupID{g},'.  '], NumSubj, 1, 700, 80);
    end
    
    parfor s = 1:length(name.SubjID{g})
        for c = 1:length(name.CondID)
            
            % Check input dataset:
            CheckInput = CheckPipelineMat(paths.InputPreprocMEG{g}{s,c}, 'ICA');
            if CheckInput == 0
                continue;
            end
            
            
            % Resample MEG data for ICA:
            cfgResampleICA           = [];
            cfgResampleICA           = FTcfg.ResampleICA;
            cfgResampleICA.inputfile = paths.InputPreprocMEG{g}{s,c};  % Read input data
            
            ResampledData = ft_resampledata(cfgResampleICA);
            
			
            % Run ICA decomposition:
            cfgICA = [];
            cfgICA = FTcfg.ICA;
			
			if exist(paths.ICAunmixing{g}{s,c}, 'file') && strcmp(RunChoice, 'APPLY')
				UnmixingInfo = LoadFTmat(paths.ICAunmixing{g}{s,c}, 'ICA');
				if isempty(UnmixingInfo)
					continue;
				end
				
				cfgICA.unmixing  = UnmixingInfo.unmixing;
				cfgICA.topolabel = UnmixingInfo.topolabel;
				
			elseif strcmp(RunChoice, 'REDO')
				if isfield(cfgICA, 'unmixing')
					cfgICA = rmfield(cfgICA, 'unmixing');
				end
				if isfield(cfgICA, 'topolabel')
					cfgICA = rmfield(cfgICA, 'topolabel');
				end
			end
			
            ICAcomponents = ft_componentanalysis(cfgICA, ResampledData);

			
            % Save ICA components, unmixing info, and call-info:
			CheckSavePath(paths.ICAunmixing{g}{s,c},   'ICA');
            CheckSavePath(paths.ICAcomponents{g}{s,c}, 'ICA');
			
            ParforSaveICA...
				(paths.ICAcomponents{g}{s,c}, paths.ICAunmixing{g}{s,c}, ICAcomponents, cfgICA);
            
            ResampledData = [];  % Free memory immediately
            ICAcomponents = [];
            
        end  % Cond
        
        if UseProgBar == 1 && mod(s, 1) == 0
            ppm.increment();  % move up progress bar
        end
        
    end  % Subj
    if UseProgBar == 1
        ppm.delete();
    end
    
end  % Group



%=========================%
% CHECK FOR OUTPUT FILES: %
%=========================%

for g = 1:length(name.GroupID)
    for s = 1:length(name.SubjID{g})
        for c = 1:length(name.CondID)
			
			if ~exist(paths.ICAunmixing{g}{s,c}, 'file')
				fprintf(ErrLog, ['ERROR: Missing ICA unmixing info file for:'...
					'\n %s \n\n'], paths.ICAunmixing{g}{s,c});
			end
            
            if ~exist(paths.ICAcomponents{g}{s,c}, 'file')
                fprintf(ErrLog, ['ERROR: Missing ICA component results file for:'...
                    '\n %s \n\n'], paths.ICAcomponents{g}{s,c});  
            end
            
        end
    end
end



%=================%

if exist([pwd,'/ErrorLog_ICA.txt'], 'file')
    LogCheck = dir('ErrorLog_ICA.txt');
    if LogCheck.bytes ~= 0  % File not empty
        open('ErrorLog_ICA.txt');
    else
        delete('ErrorLog_ICA.txt');
    end
end

fclose(ErrLog);
diary off

end



function ParforSaveICA(OutpathComp, OutpathUnmix, ICAcomponents, cfgICA)
	
	% Save ICA unmixing info separately as well:
	ICAunmixing.ICAmethod = ICAcomponents.cfg.method;
	ICAunmixing.fsample   = ICAcomponents.fsample;
	ICAunmixing.unmixing  = ICAcomponents.unmixing;
	ICAunmixing.topolabel = ICAcomponents.topolabel;
	ICAunmixing.label     = ICAcomponents.label;
	ICAunmixing.cfgICA    = cfgICA;
	
	ICAcomponents
	ICAunmixing
	
	save(OutpathUnmix, 'ICAunmixing');
    save(OutpathComp,  'ICAcomponents');
end